package fem2;

import inf.jlinalg.Array1DMatrix;
import inf.jlinalg.Array2DMatrix;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import math2.LinAlgUtilities;

/**
 * Enhanced strain mode based on classical formulation of Simo-Rifai &
 * extensions.
 * 
 * @author hbui
 * 
 */
public abstract class EASMode {

	/**
	 * @return number of enhanced strain modes
	 */
	public abstract int countEnhancedStrainModes();

	/**
	 * @return number of enhanced stress modes
	 */
	public abstract int countRecoverStressModes();

	/**
	 * compute the enhanced assumed strain interpolation matrix
	 * 
	 * @param g
	 * @param xi
	 */
	public abstract void computeEASInterpolationMatrix(IMatrix g, double... xi);

	/**
	 * compute the interpolation matrix for recover stress
	 * 
	 * @param g
	 * @param xi
	 */
	public abstract void computeRecoverStressInterpolationMatrix(IMatrix g, double... xi);

	// /**
	// * @param J
	// * @return
	// */
	// private double[] computeTransformationMatrix(int dim, double[][] J) {
	// double[] F0;
	// if (dim == 2) {
	// // F0 = new double[] { //
	// // J[0][0] * J[0][0], J[1][0] * J[0][1], 2 * J[0][0] * J[0][1], //
	// // J[0][1] * J[1][0], J[1][1] * J[1][1], 2 * J[1][0] * J[1][1], //
	// // J[0][0] * J[1][0], J[0][1] * J[1][1], (J[0][0] * J[1][1] +
	// // J[0][1] * J[1][0]) //
	//
	// F0 = new double[] { //
	// J[0][0] * J[0][0], J[0][1] * J[0][1], 2 * J[0][0] * J[0][1], //
	// J[1][0] * J[1][0], J[1][1] * J[1][1], 2 * J[1][0] * J[1][1], //
	// J[0][0] * J[1][0], J[0][1] * J[1][1], J[0][0] * J[1][1] + J[0][1] *
	// J[1][0] //
	// };
	// } else if (dim == 3) {
	// F0 = new double[] {
	// J[0][0] * J[0][0],
	// J[0][1] * J[0][1],
	// J[0][2] * J[0][2],
	// 2 * J[0][0] * J[0][1],
	// 2 * J[0][1] * J[0][2],
	// 2 * J[0][0] * J[0][2],//
	// J[1][0] * J[1][0],
	// J[1][1] * J[1][1],
	// J[1][2] * J[1][2],
	// 2 * J[1][0] * J[1][1],
	// 2 * J[1][1] * J[1][2],
	// 2 * J[1][0] * J[1][2],//
	// J[2][0] * J[2][0],
	// J[2][1] * J[2][1],
	// J[2][2] * J[2][2],
	// 2 * J[2][0] * J[2][1],
	// 2 * J[2][1] * J[2][2],
	// 2 * J[2][0] * J[2][2],//
	// J[0][0] * J[1][0],
	// J[0][1] * J[1][1],
	// J[0][2] * J[1][2],
	// J[0][1] * J[1][0] + J[0][0] * J[1][1],
	// J[0][2] * J[1][1] + J[0][1] * J[1][2],
	// J[0][2] * J[1][0] + J[0][0] * J[1][2],//
	// J[1][0] * J[2][0], J[1][1] * J[2][1],
	// J[1][2] * J[2][2],
	// J[1][1] * J[2][0] + J[1][0] * J[2][1],
	// J[1][2] * J[2][1] + J[1][1] * J[2][2],
	// J[1][2] * J[2][0] + J[1][0] * J[2][2],//
	// J[0][0] * J[2][0], J[0][1] * J[2][1], J[0][2] * J[2][2],
	// J[0][1] * J[2][0] + J[0][0] * J[2][1], J[0][2] * J[2][1] + J[0][1] *
	// J[2][2],
	// J[0][2] * J[2][0] + J[0][0] * J[2][2] //
	// };
	// } else {
	// throw new Error("invalid dimension");
	// }
	// return F0;
	// }

	/**
	 * 
	 * @param dim
	 * @param F0
	 * @param J
	 */
	private void computeTransformationMatrix(int dim, IMatrix F0, double[][] J) {
		if (dim == 2) {
			F0.set(0, 0, J[0][0] * J[0][0]);
			F0.set(0, 1, J[0][1] * J[0][1]);
			F0.set(0, 2, 2 * J[0][0] * J[0][1]);
			F0.set(1, 0, J[1][0] * J[1][0]);
			F0.set(1, 1, J[1][1] * J[1][1]);
			F0.set(1, 2, 2 * J[1][0] * J[1][1]);
			F0.set(2, 0, J[0][0] * J[1][0]);
			F0.set(2, 1, J[0][1] * J[1][1]);
			F0.set(2, 2, J[0][0] * J[1][1] + J[0][1] * J[1][0]);
		} else if (dim == 3) {
			F0.set(0, 0, J[0][0] * J[0][0]);
			F0.set(0, 1, J[0][1] * J[0][1]);
			F0.set(0, 2, J[0][2] * J[0][2]);
			F0.set(0, 3, 2 * J[0][0] * J[0][1]);
			F0.set(0, 4, 2 * J[0][1] * J[0][2]);
			F0.set(0, 5, 2 * J[0][0] * J[0][2]);
			F0.set(1, 0, J[1][0] * J[1][0]);
			F0.set(1, 1, J[1][1] * J[1][1]);
			F0.set(1, 2, J[1][2] * J[1][2]);
			F0.set(1, 3, 2 * J[1][0] * J[1][1]);
			F0.set(1, 4, 2 * J[1][1] * J[1][2]);
			F0.set(1, 5, 2 * J[1][0] * J[1][2]);
			F0.set(2, 0, J[2][0] * J[2][0]);
			F0.set(2, 1, J[2][1] * J[2][1]);
			F0.set(2, 2, J[2][2] * J[2][2]);
			F0.set(2, 3, 2 * J[2][0] * J[2][1]);
			F0.set(2, 4, 2 * J[2][1] * J[2][2]);
			F0.set(2, 5, 2 * J[2][0] * J[2][2]);
			F0.set(3, 0, J[0][0] * J[1][0]);
			F0.set(3, 1, J[0][1] * J[1][1]);
			F0.set(3, 2, J[0][2] * J[1][2]);
			F0.set(3, 3, J[0][1] * J[1][0] + J[0][0] * J[1][1]);
			F0.set(3, 4, J[0][2] * J[1][1] + J[0][1] * J[1][2]);
			F0.set(3, 5, J[0][2] * J[1][0] + J[0][0] * J[1][2]);
			F0.set(4, 0, J[1][0] * J[2][0]);
			F0.set(4, 1, J[1][1] * J[2][1]);
			F0.set(4, 2, J[1][2] * J[2][2]);
			F0.set(4, 3, J[1][1] * J[2][0] + J[1][0] * J[2][1]);
			F0.set(4, 4, J[1][2] * J[2][1] + J[1][1] * J[2][2]);
			F0.set(4, 5, J[1][2] * J[2][0] + J[1][0] * J[2][2]);
			F0.set(5, 0, J[0][0] * J[2][0]);
			F0.set(5, 1, J[0][1] * J[2][1]);
			F0.set(5, 2, J[0][2] * J[2][2]);
			F0.set(5, 3, J[0][1] * J[2][0] + J[0][0] * J[2][1]);
			F0.set(5, 4, J[0][2] * J[2][1] + J[0][1] * J[2][2]);
			F0.set(5, 5, J[0][2] * J[2][0] + J[0][0] * J[2][2]);
		}
	}

	/**
	 * compute enhanced strain operator
	 * 
	 * @param mp
	 * @param g
	 * @param xi
	 * @return
	 */
	public void computeEnhancedStrainOperator(MeshPart mp, IMatrix g, double... xi) {
		int dim = mp.getDimension();
		int ntrs = dim * (dim + 1) / 2;
		double[] center = new double[dim];
		double[][] J = mp.jacXAt(center);
		double j0 = LinAlgUtilities.determinantNxN(J);

		double[][] jac = mp.jacXAt(xi);
		double j = LinAlgUtilities.determinantNxN(jac);

		BLAM.zero(g);
		computeEASInterpolationMatrix(g, xi);
		BLAM.scale(j0 / j, g);

		// double[] F0 = computeTransformationMatrix(dim, J);

		Array1DMatrix F0 = Array1DMatrix.createRowMajor(ntrs, ntrs);
		computeTransformationMatrix(dim, F0, J);
		LinAlgUtilities.backSolve(BLAM.TRANSPOSE, F0.getValues(), g);
	}

	/**
	 * compute enhanced stress operator.<br>
	 * Note: for incompatible mode element. The enhanced stress operator is used
	 * for post process only. The stress used in computation of ri is
	 * algorithmic stress
	 * 
	 * 
	 * @param mp
	 * @param g
	 * @param xi
	 * @return
	 */
	public void computeEnhancedStressOperator(MeshPart mp, IMatrix g, double... xi) {
		int dim = mp.getDimension();
		int ntrs = dim * (dim + 1) / 2;
		double[] center = new double[dim];
		double[][] J = mp.jacXAt(center);

		Array1DMatrix F0 = Array1DMatrix.createRowMajor(ntrs, ntrs);
		computeTransformationMatrix(dim, F0, J);

		IMatrix Sigma = new Array2DMatrix(ntrs, countRecoverStressModes());
		computeRecoverStressInterpolationMatrix(Sigma, xi);

		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, F0, BLAM.NO_TRANSPOSE, Sigma, 0.0, g);

		// BLAM.zero(g);
		// double[] F0 = computeTransformationMatrix(dim, J);
		// if (dim == 2) {
		// for (int i = 0; i < 3; i++) {
		// for (int j = 0; j < 3; j++) {
		// g.set(i, j, F0[i * 3 + j]);
		// }
		// }
		// g.set(0, 3, F0[0] * xi[1]);
		// g.set(0, 4, F0[1] * xi[0]);
		// g.set(1, 3, F0[3] * xi[1]);
		// g.set(1, 4, F0[4] * xi[0]);
		// g.set(2, 3, F0[6] * xi[1]);
		// g.set(2, 4, F0[7] * xi[0]);
		// } else if (dim == 3) {
		// throw new NotImplementedError();
		// } else {
		// throw new Error("invalid dimension");
		// }
	}
}
